home *** CD-ROM | disk | FTP | other *** search
/ IRIX Base Documentation 2002 November / SGI IRIX Base Documentation 2002 November.iso / usr / share / catman / p_man / cat3 / SCSL / dlatrs.z / dlatrs
Encoding:
Text File  |  2002-10-03  |  8.9 KB  |  265 lines
  1.  
  2.  
  3.  
  4. DDDDLLLLAAAATTTTRRRRSSSS((((3333SSSS))))                                                          DDDDLLLLAAAATTTTRRRRSSSS((((3333SSSS))))
  5.  
  6.  
  7.  
  8. NNNNAAAAMMMMEEEE
  9.      DLATRS - solve one of the triangular systems  A *x = s*b or A'*x = s*b
  10.      with scaling to prevent overflow
  11.  
  12. SSSSYYYYNNNNOOOOPPPPSSSSIIIISSSS
  13.      SUBROUTINE DLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM,
  14.                         INFO )
  15.  
  16.          CHARACTER      DIAG, NORMIN, TRANS, UPLO
  17.  
  18.          INTEGER        INFO, LDA, N
  19.  
  20.          DOUBLE         PRECISION SCALE
  21.  
  22.          DOUBLE         PRECISION A( LDA, * ), CNORM( * ), X( * )
  23.  
  24. IIIIMMMMPPPPLLLLEEEEMMMMEEEENNNNTTTTAAAATTTTIIIIOOOONNNN
  25.      These routines are part of the SCSL Scientific Library and can be loaded
  26.      using either the -lscs or the -lscs_mp option.  The -lscs_mp option
  27.      directs the linker to use the multi-processor version of the library.
  28.  
  29.      When linking to SCSL with -lscs or -lscs_mp, the default integer size is
  30.      4 bytes (32 bits). Another version of SCSL is available in which integers
  31.      are 8 bytes (64 bits).  This version allows the user access to larger
  32.      memory sizes and helps when porting legacy Cray codes.  It can be loaded
  33.      by using the -lscs_i8 option or the -lscs_i8_mp option. A program may use
  34.      only one of the two versions; 4-byte integer and 8-byte integer library
  35.      calls cannot be mixed.
  36.  
  37. PPPPUUUURRRRPPPPOOOOSSSSEEEE
  38.      DLATRS solves one of the triangular systems A *x = s*b or A'*x = s*b with
  39.      scaling to prevent overflow. Here A is an upper or lower triangular
  40.      matrix, A' denotes the transpose of A, x and b are n-element vectors, and
  41.      s is a scaling factor, usually less than or equal to 1, chosen so that
  42.      the components of x will be less than the overflow threshold.  If the
  43.      unscaled problem will not cause overflow, the Level 2 BLAS routine DTRSV
  44.      is called.  If the matrix A is singular (A(j,j) = 0 for some j), then s
  45.      is set to 0 and a non-trivial solution to A*x = 0 is returned.
  46.  
  47.  
  48. AAAARRRRGGGGUUUUMMMMEEEENNNNTTTTSSSS
  49.      UPLO    (input) CHARACTER*1
  50.              Specifies whether the matrix A is upper or lower triangular.  =
  51.              'U':  Upper triangular
  52.              = 'L':  Lower triangular
  53.  
  54.      TRANS   (input) CHARACTER*1
  55.              Specifies the operation applied to A.  = 'N':  Solve A * x = s*b
  56.              (No transpose)
  57.              = 'T':  Solve A'* x = s*b  (Transpose)
  58.              = 'C':  Solve A'* x = s*b  (Conjugate transpose = Transpose)
  59.  
  60.  
  61.  
  62.  
  63.                                                                         PPPPaaaaggggeeee 1111
  64.  
  65.  
  66.  
  67.  
  68.  
  69.  
  70. DDDDLLLLAAAATTTTRRRRSSSS((((3333SSSS))))                                                          DDDDLLLLAAAATTTTRRRRSSSS((((3333SSSS))))
  71.  
  72.  
  73.  
  74.      DIAG    (input) CHARACTER*1
  75.              Specifies whether or not the matrix A is unit triangular.  = 'N':
  76.              Non-unit triangular
  77.              = 'U':  Unit triangular
  78.  
  79.      NORMIN  (input) CHARACTER*1
  80.              Specifies whether CNORM has been set or not.  = 'Y':  CNORM
  81.              contains the column norms on entry
  82.              = 'N':  CNORM is not set on entry.  On exit, the norms will be
  83.              computed and stored in CNORM.
  84.  
  85.      N       (input) INTEGER
  86.              The order of the matrix A.  N >= 0.
  87.  
  88.      A       (input) DOUBLE PRECISION array, dimension (LDA,N)
  89.              The triangular matrix A.  If UPLO = 'U', the leading n by n upper
  90.              triangular part of the array A contains the upper triangular
  91.              matrix, and the strictly lower triangular part of A is not
  92.              referenced.  If UPLO = 'L', the leading n by n lower triangular
  93.              part of the array A contains the lower triangular matrix, and the
  94.              strictly upper triangular part of A is not referenced.  If DIAG =
  95.              'U', the diagonal elements of A are also not referenced and are
  96.              assumed to be 1.
  97.  
  98.      LDA     (input) INTEGER
  99.              The leading dimension of the array A.  LDA >= max (1,N).
  100.  
  101.      X       (input/output) DOUBLE PRECISION array, dimension (N)
  102.              On entry, the right hand side b of the triangular system.  On
  103.              exit, X is overwritten by the solution vector x.
  104.  
  105.      SCALE   (output) DOUBLE PRECISION
  106.              The scaling factor s for the triangular system A * x = s*b  or
  107.              A'* x = s*b.  If SCALE = 0, the matrix A is singular or badly
  108.              scaled, and the vector x is an exact or approximate solution to
  109.              A*x = 0.
  110.  
  111.      CNORM   (input or output) DOUBLE PRECISION array, dimension (N)
  112.  
  113.              If NORMIN = 'Y', CNORM is an input argument and CNORM(j) contains
  114.              the norm of the off-diagonal part of the j-th column of A.  If
  115.              TRANS = 'N', CNORM(j) must be greater than or equal to the
  116.              infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) must be
  117.              greater than or equal to the 1-norm.
  118.  
  119.              If NORMIN = 'N', CNORM is an output argument and CNORM(j) returns
  120.              the 1-norm of the offdiagonal part of the j-th column of A.
  121.  
  122.      INFO    (output) INTEGER
  123.              = 0:  successful exit
  124.              < 0:  if INFO = -k, the k-th argument had an illegal value
  125.  
  126.  
  127.  
  128.  
  129.                                                                         PPPPaaaaggggeeee 2222
  130.  
  131.  
  132.  
  133.  
  134.  
  135.  
  136. DDDDLLLLAAAATTTTRRRRSSSS((((3333SSSS))))                                                          DDDDLLLLAAAATTTTRRRRSSSS((((3333SSSS))))
  137.  
  138.  
  139.  
  140. FFFFUUUURRRRTTTTHHHHEEEERRRR DDDDEEEETTTTAAAAIIIILLLLSSSS
  141.      A rough bound on x is computed; if that is less than overflow, DTRSV is
  142.      called, otherwise, specific code is used which checks for possible
  143.      overflow or divide-by-zero at every operation.
  144.  
  145.      A columnwise scheme is used for solving A*x = b.  The basic algorithm if
  146.      A is lower triangular is
  147.  
  148.           x[1:n] := b[1:n]
  149.           for j = 1, ..., n
  150.                x(j) := x(j) / A(j,j)
  151.                x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]
  152.           end
  153.  
  154.      Define bounds on the components of x after j iterations of the loop:
  155.         M(j) = bound on x[1:j]
  156.         G(j) = bound on x[j+1:n]
  157.      Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.
  158.  
  159.      Then for iteration j+1 we have
  160.         M(j+1) <= G(j) / | A(j+1,j+1) |
  161.         G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |
  162.                <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )
  163.  
  164.      where CNORM(j+1) is greater than or equal to the infinity-norm of column
  165.      j+1 of A, not counting the diagonal.  Hence
  166.  
  167.         G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )
  168.                      1<=i<=j
  169.      and
  170.  
  171.         |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )
  172.                                       1<=i< j
  173.  
  174.      Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTRSV if the
  175.      reciprocal of the largest M(j), j=1,..,n, is larger than
  176.      max(underflow, 1/overflow).
  177.  
  178.      The bound on x(j) is also used to determine when a step in the columnwise
  179.      method can be performed without fear of overflow.  If the computed bound
  180.      is greater than a large constant, x is scaled to prevent overflow, but if
  181.      the bound overflows, x is set to 0, x(j) to 1, and scale to 0, and a
  182.      non-trivial solution to A*x = 0 is found.
  183.  
  184.      Similarly, a row-wise scheme is used to solve A'*x = b.  The basic
  185.      algorithm for A upper triangular is
  186.  
  187.           for j = 1, ..., n
  188.                x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)
  189.           end
  190.  
  191.      We simultaneously compute two bounds
  192.  
  193.  
  194.  
  195.                                                                         PPPPaaaaggggeeee 3333
  196.  
  197.  
  198.  
  199.  
  200.  
  201.  
  202. DDDDLLLLAAAATTTTRRRRSSSS((((3333SSSS))))                                                          DDDDLLLLAAAATTTTRRRRSSSS((((3333SSSS))))
  203.  
  204.  
  205.  
  206.           G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j
  207.           M(j) = bound on x(i), 1<=i<=j
  208.  
  209.      The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we add
  210.      the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.  Then the
  211.      bound on x(j) is
  212.  
  213.           M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |
  214.  
  215.                <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )
  216.                          1<=i<=j
  217.  
  218.      and we can safely call DTRSV if 1/M(n) and 1/G(n) are both greater than
  219.      max(underflow, 1/overflow).
  220.  
  221.  
  222. SSSSEEEEEEEE AAAALLLLSSSSOOOO
  223.      INTRO_LAPACK(3S), INTRO_SCSL(3S)
  224.  
  225.      This man page is available only online.
  226.  
  227.  
  228.  
  229.  
  230.  
  231.  
  232.  
  233.  
  234.  
  235.  
  236.  
  237.  
  238.  
  239.  
  240.  
  241.  
  242.  
  243.  
  244.  
  245.  
  246.  
  247.  
  248.  
  249.  
  250.  
  251.  
  252.  
  253.  
  254.  
  255.  
  256.  
  257.  
  258.  
  259.  
  260.  
  261.                                                                         PPPPaaaaggggeeee 4444
  262.  
  263.  
  264.  
  265.